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ABSTRACT 

Full  field  numerical  solutions  for  a  crack  which  lies  along  the  interface  of  an 
elastic-plastic  medium  and  a  rigid  substrate  are  presented.  The  solutions  are 
obtained  using  a  small  strain  version  of  the  J2  deformation  theory  with  power  law 
strain  hardening.  In  the  present  article  results  for  loading  causing  only  small  scale 
yielding  at  the  crack  tip  are  described;  in  a  subsequent  article  results  for  contained 
yielding  and  fully  plastic  behavior  will  be  presented.  The  oscillatory  stresses  on 
the  bond  line  and  overlapping  of  the  crack  faces,  characteristic  of  small  strain  linear 
elasticity,  are  shown  to  be  essentially  precluded  by  material  nonlinearity.  In 
addition  we  find  that,  although  the  near  tip  fields  do  not  appear  to  have  a  separable 
form  as  for  the  well  known  HRR  fields  in  homogeneous  media,  they  do  bear 
interesting  similarities  to  certain  mixed  mode  HRR  fields.  Numerical  procedures 
appropriate  for  solving  a  general  class  of  interface  crack  problems  are  described. 
A  description  of  a  numerical  method  for  extracting  the  mixed  mode  stress  intensities 
for  cracks  at  interfaces,  and  in  homogeneous  isotropic  or  anisotropic  media,  is 


included. 


1.  INTRODUCTION 


Fracture,  whether  it  occurs  by  fibrous  or  cleavage  modes,  invariably  begins  at 
interfaces.  In  structural  alloys,  precipitation  or  segregation  of  impurities  at  grain  boundaries 
can  lead  to  transitions  from  ductile  to  cleavage  fracture  and  to  concomitant  losses  in 
ductility  and  toughness.  In  so-called  advanced  materials  such  as  structural  ceramics, 
composites,  and  polycrystalline  intermetallic  alloys,  interfacial  and  intergranular  fractures 
are  common  and  may,  in  large  part,  determine  the  material’s  overall  mechanical  response.  At 
present,  however,  what  exists  is  only  a  relatively  small  group  of  solutions  for  the  crack  tip 
fields  at  the  interfaces  of  isotropic  linear  elastic  media  (e.g.  Williams,  1959;  Erdogan,  1963, 
1965,  Sih  and  Rice,  1964;  Rice  and  Sih,  1965;  England,  1965).  A  formal  solution  for  the 
field  of  an  isolated  crack  lying  along  the  interface  of  linear  anisotropic  media  has  been  given 
by  Willis  (1971)  which,  while  perhaps  not  lending  itself  to  a  ready  evaluation  of  the  nature  of 
the  field,  does  allow  for  a  computation  of  the  energy  release  rate  for  crack  extension  along 
the  interface.  To  our  knowledge  there  are  no  solutions  for  the  fields  of  interface  cracks 
in  elastic-plastic  materials  and  this  makes  it  difficult  to  analyze  the  legion  of  interface 
fracture  phenomena  alluded  to  above.  One  purpose  of  this  study  is  to  provide  such  solutions 
within  a  framework  consistent  with  what  exists  for  nonlinear  fracture  mechanics  in 
homogeneous  materials  (Hutchinson,  1983).  Another  is  to  develop  numerical  procedures 
capable  of  accurately  resolving  the  complex  fields  which  develop  at  interfaces.  Our 
present  >n  is  to  be  given  in  two  parts:  in  the  present  paper  the  structure  of  the  asymptotic 
field  under  small  scale  yielding  conditions  is  described;  the  fields  that  develop  under 
contained  yielding  and  fully  plastic  conditions  will  be  presented  in  a  subsequent  publication. 
Implications  of  the  characteristic  stress  and  deformation  fields  that  form  at  and  near 
interfaces  for  fracture  mechanics  will  be  detailed  in  the  latter. 

Specifically,  the  plane  strain  boundary  value  problem  for  a  center  cracked  plate  loaded 
with  uniform  remote  stresses  is  solved  and  discussed  here.  This  is  illustrated  along  with  a 
version  of  the  finite  element  mesh  used  for  numerical  solution  in  Fig.  3.  The  crack  lies 


-3- 


on  the  interface  between  a  rigid  substrate  and  a  nonlinear  elastic-plastic  material 
characterized  by  a  J2  deformation  theory.  This  problem  serves  as  a  prototype  for  a  wide 
class  of  physically  interesting  phenomena  involving  microcracking  and  macroscopic 

fracture,  while  at  the  same  time  allowing  for  an  unambiguous  and  definitive  study  of 
the  structure  of  the  fields  of  interface  cracks  within  the  context  of  a  complete  boundary 
value  problem.  Crack  tip  fields  have  been  evaluated  numerically  for  linear  elastic 
materials  and  nonlinear  elastic-plastic  materials  (characterized  by  a  small  strain  J2 
deformation  theory)  under  small  scale  yielding,  contained  yielding,  and  f ully-plastic 

conditions.  In  the  present  paper  attention  is  focused  on  clastic  behavior  and 
elastic-plastic  behavior  under  small  scale  yielding.  For  this  analysis  a  domain 
formulation  of  an  interaction  energy  release  rate  is  employed  to  extract  the  complex 
stress  intensity  factors  from  the  numerical  fields.  Some  discussion  of  results  for 
contained  yielding  is  provided,  but  as  noted  above,  a  full  treatment  of  these  cases  along 
with  the  fields  for  cracks  under  loading  conditions  that  cause  fully  yielded  behavior  will 
be  given  in  a  subsequent  article. 

In  the  context  of  interface  cracks,  small  scale  yielding  is  meant  to  pertain  to 

loading  conditions  for  which  the  plastic  zone  emanating  from  the  crack  tip  is  surrounded 
by  elastic  fields  which  are  well  approximated  by  the  elastic  singularity  fields  for  the 
interface  crack.  We  do  not  suggest  that  the  plastic  zones  develop  in  a  self-similar 

manner  as  they  do  for  cracks  in  homogeneous  media.  Indeed  the  change  in  plastic  zone 
size  and  shape  with  increasing  remote  stress  is  detailed  in  Sections  4.3,  and  this  should 
be  kept  in  mind  in  interpreting  the  results  presented  in  Sections  4.3  and  5. 

As  mentioned,  the  calculations  are  performed  for  a  material  described  by  a  small 
strain,  isotropic  J2  deformation  theory.  This  was  done  in  order  to  facilitate  the 
connection  between  the  interface  crack  solutions,  the  existing  framework  of  nonlinear 
fracture  mechanics,  and  specific  solutions  for  crack  tip  fields  in  homogeneous  media. 
Both  qualitative  and  quantitative  similarities  between  the  fields  are  highlighted  in  the 
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results  we  present.  It  is  found,  for  example,  that  the  near-tip  stress  and  strain  fields 
are  inherently  mixed  mode  regardless  of  the  far-field  loading.  Although  the  fields  for 
the  interface  crack  are  not  of  the  usual  separable  form,  there  are  nonetheless  strong 
similarities  between  them  and  the  mixed  mode  fields  for  cracks  in  homogeneous  media. 
(Thus  far,  our  attempts  to  find  fields  near  the  tip  of  an  interface  crack  which  are  of  a 
separable  form  have  not  been  successful.)  These  features  have  important  implications  for 
the  mathematical  characterization  of  the  fields  as  well  as  for  failure  processes  that  may 
be  induced  at  or  near  the  interface. 

We  note,  however,  that  the  stress  histories  that  develop  near  the  tip  of  the  interface 
crack  cannot  be  described  as  proportional  stressing.  At  low  remote  stress  levels,  for 
example,  the  fields  that  develop  near  the  tip  are  of  a  strongly  mixed  mode 
nature  -  reflecting  the  mixed  mode  character  of  the  linear  elastic  field.  At  higher 
remote  stresses  our  results  show  that  the  crack  tip  fields  are  more  nearly  like  those 
corresponding  to  pure  mode  I  fields  of  cracks  in  homogeneous  media.  The  transition 
from  strongly  mixed  mode  fields  to  nearly  mode  I  fields  certainly  involves 
non-proportional  stressing.  Thus  interface  crack  tip  behavior  is  more  complex  than  is 
found  for  cracks  in  homogeneous  media  for  which  proportional  remote  loading  induces 
proportional  stressing  near  the  crack  tip.  The  link  between  the  fields  calculated  here 
for  a  T2  deformation  theory  material  and  those  appropriate  to  path-dependent 
incremc.  al  elastic-plastic  theories  is  the  subject  of  the  current  investigation. 

The  plan  of  the  paper  is  as  follows.  In  the  next  section  we  present  a  summary  of 
available  results  for  crack  tip  stress  and  displacement  fields  in  linear  isotropic  media. 
This  introduces  the  notion  of  a  stress  concentration  vector,  or  stress  intensity  factors,  for 
interface  cracks  and  provides  the  requisite  background  to  analyze  the  elastic-plastic 
results.  A  complete  description  of  the  boundary  value  problem  to  be  solved  is  also  given. 
Numerical  methods  are  then  described,  followed  by  the  presentation  of  the  full  field 
solutions.  A  short  discussion  on  the  nature  of  the  computed  fields  concludes  this  paper. 
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2.  PERSPECTIVES  ON  CRACKS  AT  BIMATERIAL  INTERFACES 
2-1  Overview  of  Results  for  Linear  Elastic  Merlin 

A  number  of  solutions  for  the  stress  and  displacement  fields  for  cracks  lying  along 
bimaterial  interfaces  have  been  obtained  for  isotropic  materials  by,  for  example,  Williams 
(1959),  England  (1965),  Erdogan  (1963,1965),  Sih  and  Rice  (1964)  and  Rice  and  Sih  (1965). 
For  anisotropic  materials  Willis  (1971)  has  provided  a  formal  solution  which  allows  the 
energy  release  rate  to  be  evaluated.  Our  concern  in  this  section  is  to  illustrate  the 
nature  of  the  tractions  and  the  displacement  discontinuities  on  the  crack  line  at  the 
crack  tip,  identify  stress  intensity  factors  and  relate  them  to  energy  release  rates.  The 

discussion  is  confined  to  isotropic  elastic  media  where  the  results  are  easily  represented 
in  a  concise  structure. 


We  refer  to  Fig.  1  which  shows  a  crack  lying  along  an  interface  separating 
media  1  and  2.  The  shear  moduli  and  Poisson’s  ratios  are  /zl5  /z2  and  Vp  v2.  ?  is  the 

distance  measured  along  the  bond  line  from  the  crack  tip  such  that  £  =  r,  the  polar 

coordinate,  when  9  =  0  and  Q  =  n.  In  what  follows  it  is  to  be  understood  that  ?  is 

directed  along  the  bond  line  in  the  expressions  for  tractions,  whereas  it  is  directed  back 

along  the  crack  faces  in  the  expressions  given  below  for  the  displacement  jumps  across 
the  crack  face.  Existing  solutions  show  that  the  tractions  on  the  bond  line  can  be 


expressed  in  the  form 

KjCoste/nU^a))  K2sin(e/n(?/2a)) 

^22  =  —  “  • - 

✓277?  ^277? 


(2.1a) 


and 


o 


12 


K1sin(€/n(?/2a))  K2cos(e/n(?/2a)) 

- : - -  +  - - - 

✓277?  y277? 


(2.1b) 


where  e  is  the  bimaterial  constant  introduced  by  Williams  (1959),  Erdogan  (1963),  Sih  and 
Rice  (1964),  and  England  (1965),  defined  as 
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3-4v 


l  J_ 
+  ^2 
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3-4v„ 


(2.2) 


Ki  and  K2  define  two  stress  intensity  factors  which,  we  note,  have  units  of 
stress-(length)^  in  analogy  with  the  definition  of  stress  intensity  factors  in  homogeneous 
media.  It  is  convenient  to  introduce  a  complex  stress  intensity  factor  defined  by 
K  =  Kx  +  iK2.  This  definition  is  related  in  a  simple  way  to  others  available  in  the 
literature  and,  for  example,  is  related  to  the  complex  stress  intensity  factor  of  Sih  and 
Rice  (1964),  k  =  kj  +  ik2,  by  K  =  -jn  cosh(77e)(2a)*ek. 

Combining  (2.1)  as  in  the  original  solutions  cited  above  (see  also  Willis’  (1971)  eq.  (4.4)), 
the  tractions  on  the  bonding  surface  can  also  be  expressed  as  the  complex  vector 


t  °22  +  i  o12 


K 

✓2775 


ei€/n  (5/  2a)  _  _ ^  (f /?a)^ 
✓277? 


(2.3) 


Explicit  examples  of  K  for  three  crack  geometries  and  loading  configurations,  taken 
from  the  above  references,  are  given  in  Fig.  2.  We  note  that,  in  general,  the  two 
stress  modes  defined  by  Kx  and  K2  each  involve  mixed  tension  and  shear  on  the 
bonding  surface  and,  in  isotropic  media,  represent  orthogonal  modes.  That  is,  if 


then 


Ki 

—  1  eie/n(5/2a) 

✓2775 


(2.4a) 


t 


2 


K2 

✓2775 


ei{e/n  (5/2a)  +  77/2} 


(2.4b) 


which  is  oriented  at  a  right  angle  to  tx  when  plotted  on  a  graph  with  a22  on  the  real  axis 
and  o12  on  the  imaginary  axis.  When  the  two  media  are  identical,  e  =  0  and  Kx  and 

K2  reduce  to  the  usual  mode  I  and  II  stress  intensity  factors,  Kj  and  Kjj. 
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Displacement  jumps  across  the  crack  face  have  also  been  given  (explicitly  by 
England,  1965)  and,  at  the  crack  tip,  take  the  form. 


AU  =  Au2  +  i  AUj  = 


This  can  also  be  written  as 


l  H2  J 


(l+2ie)  coshne 


Kvt/2  n  eie/n^/2a) 


(2.5) 


AU  = 


1  H  [ 

(l+4e2)^  cosh7t€ 


K  v*,/2n(  5/2a)iee"i6 


(2.6) 


where  0  =  tan_12e.  Thus  the  contributions  to  AU  associated  with  the  stress  modes  AT 
and  K2  (call  them  AUj  and  AU2)  are  also  orthogonal  but  are  not  precisely  aligned  with  tj 
and  t2.  In  fact,  if  tx  and  t2  are  placed  on  a  graph  whose  axes  are  o22  and  a12  then  the 
orthogonal  pair  AU1?  AU2,  when  graphed  with  respect  to  axes  Au2  and  AUj,  would  appear 
to  be  rigidly  rotated  with  respect  to  the  orthogonal  pair  tv  t2  by  the  angle  0. 

As  noted  by  Williams  (1959),  and  in  subsequent  works,  a  characteristic  of  the  linear 
elastic  fields  of  interface  cracks  is  that  they  violate  compatibility  at  distances  very  close 
to  the  tip,  the  stresses  and  displacement  fields  on  the  crack  line  oscillate  with  unbounded 
amplitude  and  vanishing  wavelength.  Looking  ahead  to  the  problem  of  the  center  cracked 
panel  under  far  field  tension  that  we  investigate  numerically,  we  note  that  Rice  and 
Sih’s  (1965)  results  yield  the  stress  intensity  vector  K  =  o“2  vnsi  (l+2ie).  This  means, 
for  example,  that  the  normal  tractions  on  the  bond  line  become  compressive  as  soon  as 


Rejo+2ie)eie/n(S/2a)j  <  0  > 

or  using  the  definition  for  0  introduced  earlier,  when 
Re  jei[£^n  U/2a)+£]  j  q 


(2.7) 


(2.8) 
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Since  (§/2a)  <  1,  this  first  occurs  as  soon  as 


€/n(?/2a)  +  B  <  -n/2  .  (2.9) 

Thus  egg  first  becomes  compressive  at 

U/2a)  =  e“[77/2+,3]/€  .  (2.10) 

For  the  case  of  a  crack  on  a  rigid  substrate,  where  in  the  deformable  half  space  v  =  0.3, 
e  =  0.0935,  and  this  leads  to  a  value  of  5/2a  ^  7xi(r9.  Our  numerical  results 
reproduced  this  (albeit  non-physical)  feature  of  the  linear  elastic  field.  In  fact,  looking 
ahead  to  Fig.  8,  where  the  numerical  solution  for  the  interface  crack  in  a  linear  elastic 
medium  is  presented,  it  may  be  noted  that  along  a  ray  at  an  angle  of  7.5°  with  the  bond 
line  (i.e.  essentially  on  the  bond  line)  o@0  becomes  compressive  at  r/a  =  lxlO"9  which  is 
completely  consistent  with  this  stress  component  first  becoming  compressive  at 
5/a  =  14xl0'9  at  0  =  0°.  However,  we  note  that  when  the  loading  is  a  combined  remote 
tension  and  shear,  K  =  [(a22  -  2o12e)  +  i(o”2  +  2eo”2)]1/7Ta  and  a  similar  analysis  shows 
that  a00  first  becomes  compressive  on  the  bond  line  when 


(5/2a)  =  +  ^]/6 


(2.11) 


where 


X 


Now 


CO 


suppose 

00  \ 
°2  2’  ^ 


=  tan  V”  +  2eo” )/(a”  -  2ea~)}.  (2.12) 

o”2  >  0  but  o”2  =  -2o“2;  then  X  =  -58.7°  and  5/2a  ^  lxicr3  !  If 

=  -38.2  and  5/2a  =  3xl0‘5.  Thus  the  regions  where  the  stresses  (based 


on  small  strain  linear  elasticity  assumptions)  oscillate  are  not  always  confined  to  vanishingly 
small  regions  near  the  crack  tip  and  their  extents  are  quite  sensitive  to  the  nature  of  the 
remote  loading  and  the  value  of  €.  In  Section  5  full  field  results  are  presented  for  the 
interface  crack  in  nonlinear  materials  subjected  to  remote  mixed  loadings.  They  demonstrate 
that  such  oscillations  are  not  present  when  the  remote  stresses  are  small  fractions  of  the 
material  s  yield  stress.  This  is  also  shown,  in  detail,  in  Section  4  for  the  case  of  pure  remote 
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tensile  loading.  It  may  be  noted  that  when  v  =  1/2,  €  =  0  and  the  near-tip 

displacements  vanish  identically  on  the  bond  line  for  the  plane  strain  problem.  For  this 
case  the  fields  are  non-oscillatory  and  are  identical  to  the  fields  for  the  corresponding 
homogeneous  media.  This  fact  has  been  noticed  by  Knowles  and  Sternberg  (1983). 

When  written  out  in  full  the  asymptotic  crack  tip  stress  field  is  given  by 

2 

V  ( a) 

°ij  =  )  - 77  8ij  ( 9,  e/n(r/2a);e]  (2.13) 

L  (2nr)^ 
a=l 

(a) 

and,  for  completeness,  the  functions  gjj  are  given  in  Appendix  I  in  a  form  consistent 
with  our  definition  of  the  complex  stress  intensity  vector  K.  It  should  be  noted  that 
the  fields  (2.13)  are  not  of  a  separable  form.  In  contrast,  the  singular  fields  of  cracks  in 
homogeneous  media  are  separable.  It  should  also  be  noted  that  with  our  definition  of 
the  complex  stress  intensity  factor,  crack  geometry  appears  in  the  angular  functions  as 
well  as  in  the  stress  intensity  factors. 

The  energy  release  rate  is  readily  computed  from  the  given  expressions  for  the 
traction  and  displacement  fields  as  the  virtual  work  integral 

Aa 

c"Lm-0  i  I  tc«;0).aii(5;0+iaxu  (2.14) 

0 

where  t(£;0)  is  the  complex  traction  vector  at  the  crack  tip,  given  by  (2.3),  when  the 
crack  is  at  x  =  0,  and  AU(£;  0+Aa)  is  the  complex  conjugate  to  the  displacement  jump 
vector  when  the  crack  has  advanced  to  the  position  0+Aa.  The  integral  is  easily  reduced 
to  the  form 


1 

(1  -2i€) 


J*-ie 


0 


2  77  cosh  77€ 


X 


dt 


(2.15) 
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where  the  last  integral  is  recognizable  as  the  complex  Beta  function  B(l/2+ie,  3/2-ie). 
Upon  evaluation  of  B, 


G 


(2.16) 


4cosh2  77€ 


Orthogonality  of  the  stress  modes  guarantees  that  they  decouple  such  that  there  are  no 
terms  of  the  form  Note  also  that  the  integral  in  (2.14)  has  no  imaginary  part,  a 

fact  noticed  in  the  general  anisotropic  case  by  Willis  (1971). 

List  of  stress  intensity  factors 

We  complete  this  section  with  the  short  list  of  stress  intensity  factors  tabulated  in 
Fig.  2.  Included  in  this  list  is  the  center  cracked  panel  we  have  studied  numerically.  It 
may  be  noted  in  these  examples  that  with  K  defined  by  (2.3),  the  ratio  of  Kx  to  K2  does 
not  depend  on  the  characteristic  dimension  2a. 

2.2  The  Interface  Crack  Small  Scale  Yielding  Problem 

The  specific  problem  we  are  concerned  with  here,  illustrated  in  Fig.  3,  is  that  of  a 
large  plate  loaded  remotely  by  uniform  stresses.  Crack  length,  2a,  is  the  characteristic 
dimension  with  respect  to  which  the  stress  intensity  factors  listed  in  Fig.  2,  and  the 
asymptotic  angular  stress  functions  given  in  Appendix  I,  are  evaluated.  Uniform  normal 
stresses  are  imposed  on  the  external  boundaries  with  the  primary  variable  being  o”2.  No 
remote  shear  is  imposed  although,  as  is  evident  from  the  elastic  solutions,  large  shear 
stresses  develop  on  and  near  the  crack  line  due  to  the  differences  in  material  properties. 
This  is  also  true  for  the  elastic-plastic  cases  treated  numerically. 

As  noted  by  Rice  and  Sih  (1965),  continuity  of  the  extensional  strain  en  across  the 
bond  line  requires  that  the  normal  stress  parallel  to  the  bond  line  be  discontinuous. 
They  derived  the  required  jump  in  o!°.,  viz. 
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where 


co  oo  (3+w)e277e  (3w+l) 

Jn  (2)  =  wo  (1)  +  - — 

i  .  „  2  77€ 


M2(1-vx) 


(2.17) 


(2.18) 


For  the  present  case  of  the  crack  on  the  rigid  substrate,  w  -*  00  and  the  above  relation 
takes  the  form 


'll 


(1)  = 


3-e  2  77€ 


1  +  e 


2  77  e 
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(2.19) 


In  the  numerical  solutions  described  next  the  normal  stress  given  by  (2.19),  for  v  =  0.3, 
was  imposed  on  the  side  faces  of  the  deformable  medium  shown  in  Fig.  3. 

The  small  scale  yielding  study  is  carried  out  within  the  context  of  small  strain 
theory.  The  deformable  medium  is  taken  to  be  described  by  J2  deformation  theory  for  a 

Ramberg-Osgood  stress-strain  behavior.  In  uniaxial  tension  the  material  deforms 
according  to 

e/€o  =  °/ao  +  “(a/a0)n  (2.20) 


where  aQ  and  6Q  are  the  reference  stress  and  strain,  a  is  a  material  constant  (taken  to  be 
0.1  in  this  study)  and  n  is  the  strain  hardening  exponent.  Under  multi-axial  stress 


states  ay  the  strain  is 


•ij 


1+v 

- S::  + 
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l-2v 


3E  KK  2 


n— 1  s 


ij 


(2.21) 


Here,  Sjj  is  the  stress  deviator,  aQ  =  (3sjjSjj/2)^  is  the  effective  stress,  and  v  and  E  are 
the  elastic  constants.  In  writing  (2.21)  we  have  used  the  connection  aQ  =  €QE. 
Numerical  solutions  under  plane  strain  assumptions  have  been  obtained  for  a  range  of 


hardening  exponents.  However,  only  results  for  n=3  (a  high  hardening  material)  and 
n=10  (a  moderate  to  low  hardening  material)  will  be  presented.  In  all  the  calculations  we 
have  used  v  equal  to  0.3. 
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3.  NUMERICAL  PROCEDURE 

3.1  Selective  Integration 

The  standard  displacement-based  serendipity  and  Lagrangian  isoparametric 
quadrilateral  elements  perform  poorly  as  the  deformation  progresses  into  the  fully  plastic 
(incompressible)  range  (Nagtegaal  et  al.,  1974).  These  difficulties  can  be  alleviated  by 
selective/reduced  integration  and  one  implementation  of  this  is  the  so-called  B-bar 
method  (Hughes,  1980).  The  strain-displacement  matrix,  referred  to  as  the  B  matrix,  is 
assembled  in  the  following  manner.  The  deviatric  components  of  the  B  matrix  are 
evaluated  at  the  regular  quadrature  points.  Volumetric  components  are  evaluated  at  the 
reduced  quadrature  points;  the  volumetric  B  matrix  associated  with  the  regular 
quadrature  points  is  obtained  by  interpolation/extrapolation  using  a  lower  order  shape 
function.  The  volumetric  and  deviatoric  components  are  then  combined  to  give  the 
desired  matrix  associated  with  the  regular  quadrature  points.  The  latter  matrix,  denoted 
by  B  (to  distinguish  it  from  the  usual  B  matrix,  Zienkiewicz,  1977),  is  used  in  the 
formation  of  element  stiffness  matrix  and  for  evaluating  strains. 

The  full  Lagrangian  shape  functions  are  necessary  for  approaching  the 
incompressible  limit  (Malkus  and  Hughes,  1978).  In  this  regard  it  appears  that  an 
optimal  element  for  treating  nearly  incompressible  deformation  is  the  9-node  Lagrangian 
quadrilateral  element.  For  this  element  the  deviatoric  part  of  B  is  evaluated  at  3x3 
Gaussian  quadrature  points  while  the  volumetric  part  is  evaluated  at  2x2  quadrature 
points.  The  volumetric  B  at  each  of  the  3x3  quadrature  points  is  obtained  by  bilinear 
interpolation/extrapolation  and  combined  with  the  deviatoric  part  to  give  B,  i.e. 

B  =  Bdev  +  Bvo1  .  (3  n 

The  (tangent)  stiffness  matrix  for  an  element  is  given  by  the  sum  of  the  inner  products 
at  each  Gauss  point. 


where  D  is  the  matrix  of  material  moduli  evaluated  at  the  3x3  quadrature  points  and  Wp 
is  the  appropriate  weight.  It  is  known  that  the  selective  integration  method  gives  rise  to 
a  spurious  pressure  mode  in  the  4-node  quadrilateral  element.  For  this  element,  Bvo1  is 
evaluated  with  one-point  quadrature.  Belytschko  and  Tsay  (1983)  have  proposed  a 
stabilization  procedure  to  suppress  such  spurious  modes.  Spurious  pressure  modes  are 
avoided  by  using  the  9-node  Lagrangian  element. 

We  obtained  results  that  are  nearly  identical  in  our  numerical  experimentations 
with  the  B-bar  method  and  selective  integration  (2x2  quadrature  on  volumetric  stiffness 
matrix  and  3x3  quadrature  on  deviatoric  stiffness  matrix)  on  several  boundary  value 
problems  based  on  the  volume  preserving  plasticity  relations  (2.20).  In  this  connection,  it 
may  be  noted  that  accurate  solutions  to  fully  plastic  crack  problems  in  homogeneous 
media  have  been  obtained  using  the  9-node  Lagrangian  element  in  conjunction  with 
selective  integration  and  parameter  tracking  (e.g.  Shih  and  Needleman,  1984).  For  the 
present  boundary  value  problem  based  on  (2.21),  the  incompressibility  constraint  is  more 
directly  and  efficiently  accommodated  by  the  B-bar  method.  In  passing,  we  point  out 
that  the  results  plotted  in  the  figures  of  this  paper  are  the  actual  values  computed  at  the 
quadrature  points  -  smoothing  has  not  been  applied  to  any  data. 

3.2  Parameter  Tracking 

The  solution  to  the  nonlinear  boundary  value  problem  is  obtained  by  the 
Newton-Raphson  method.  The  iterative  method  is  second  order  convergent  if  a  close 
initial  estimate  of  the  solution  is  available.  For  highly  nonlinear  problems  of  the  type 
being  investigated  here,  the  initial  estimate  is  generated  by  parameter  tracking.  We  begin 
by  obtaining  the  solution  to  the  linear  elastic  problem  at  the  desired  remote  stress.  The 
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linear  elastic  solution  is  then  employed  as  the  initial  estimate  in  the  iterations  for  the 
mildly  nonlinear  problem  with  n=2.  The  convergent  solution  for  the  mildly  nonlinear 
problem  is  then  employed  as  the  initial  estimate  for  a  more  nonlinear  problem,  say  n=3. 
In  this  manner  solutions  for  material  characteristics  ranging  from  high  hardening  to 
nearly  non-hardening  are  obtained.  Similar  tracking  is  performed  on  the  remote  stress 
levels;  that  is  once  a  convergent  solution  is  obtained  for  a  given  remote  loading  and 
nonlinear  material  description,  it  may  be  used  as  the  initial  guess  for  a  new  solution  at  a 
higher  remote  load. 

3.3  Domain  Representation  of  J 

We  consider  the  line  integral  (Rice,  1968)  defined  by 


(3.3) 


Here,  r  is  a  contour  beginning  at  the  bottom  crack  face  and  ending  on  the  top  face  and 
nj  is  the  outward  normal  to  I"  as  shown  in  Fig.  4.  For  a  nonlinear  elastic  solid  or 
deformation  theory  solid,  the  integrand  is  divergence  free.  With  regard  to  the  crack 
between  a  deformable  medium  and  a  rigid  substrate,  the  contour  r  begins  on  the  bonding 
plane  and  ends  on  the  upper  crack  face.  Since  there  are  no  contributions  to  the  integral 
along  the  traction  free  crack  face  and  the  interface  Xj  =  0"*^  (uj  vanishes  on  the  bond 
line),  t  le  value  of  J  does  not  depend  on  the  contour  r,  i.e.  the  integral  is 
path-independent.  For  the  general  bimaterial  problem,  the  contour  r  begins  on  the 
lower  crack  face  and  ends  on  the  upper  crack  face.  Path-independence  is  easily 
demonstrated  for  this  general  case  once  it  is  recognized  that  the  contribution  from  the 
upper  interface,  X2  =  0  ,  negates  the  contribution  from  the  lower  interface,  X2=  0  For 
monotonic  loading  conditions  and  proportional  stress  history  at  each  material  point,  the 
above  integral  is  also  path-independent  for  elastic-plastic  solids. 


For  numerical  purposes  it  is  more  advantageous  to  recast  the  line-integral  (3.3)  as 
an  area/domain  integral.  Such  a  representation  is  naturally  compatible  with  the  finite 
element  method  and  very  accurate  values  for  J  have  been  obtained  using  the  domain 
representation.  A  domain  integral  formulation  for  the  crack  tip  force  (including  the 
J-integral)  for  general  material  response  and  arbitrary  crack  tip  motion  has  been  detailed 
by  Moran  and  Shih  (1986)  and  has  been  applied  to  several  crack  problems  in 
homogeneous  media  (e.g.  Li  et  al.,  1985,  Shih  et  al.,  1986).  To  obtain  the  desired 
domain  representation  for  J,  weighting  functions  qj  are  introduced.  For  two 
dimensional  problems  with  the  crack  line  oriented  in  the  x1-direction,  qj  is  the  only 
non-zero  function  and  it  has  the  value  of  unity  on  the  contour  r  and  zero  on  the  outer 
contour  CQ  shown  in  Fig.  4.  Within  the  area  enclosed  by  r,  CQ  and  the  crack  faces  C+ 
and  C_,  qj  is  an  arbitrary  smooth  function  of  Xj^  and  x2  with  values  in  the  range  from 
zero  to  one.  The  function  q1(x1,x2)  may  be  interpreted  as  the  virtual  translation  of  the 
material  point  at  (xrx2)  due  to  a  unit  virtual  extension  of  the  crack  in  the 
x1-direction.  Development  along  these  latter  lines  has  been  given  by  Parks,  1977,  Hellen, 
1975  and  deLorenzi,  1982. 

For  a  function  qx  that  satisfies  the  above  conditions,  the  integral  in  (3.3)  can  be 
restated  as  a  line  integral  over  the  closed  contour  C  which  consists  of  F,  CQ  and  the 
crack  faces  defined  by 


J  = 

• 

3uj 

aT  -  W6ii 

c 

k.  L  „ 

m 


j  QjdC. 


(3.4) 


Here  rnj  are  the  components  of  the  unit  vector  normal  to  C  that  points  away  from  the 
enclosed  area.  On  the  contour  r,  mj  is  equal  to  the  negative  of  nj  which  has  been 
defined  earlier.  The  crack  faces  are  assumed  to  be  traction  free.  Applying  the 
divergence  theorem  to  (3.4)  and  using  the  equation  of  equilibrium,  the  equivalent  domain 
representation  of  (3.3)  is 


(3.5) 


J  = 


9ui 

**  iT  "W8u- 


9Qj 
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where  A  is  the  area  enclosed  by  the  closed  contour  C. 
is  taken  to  be  a  potential  function  of  stress,  i.e.,  a-j  = 


In  deriving  the  expression  (3.5),  W 
9W/06JJ. 


3.4  Finite  Element  Representation  of  Domain  J 

The  finite  element  form  of  the  domain  representation  (3.5)  is  briefly  discussed.  In 
isoparametric  (771,r?2)  (-1  $  *?•  ^  1,  1=1,2)  space,  the  shape  functions  are  constructed  from 
the  basic  functions 

1  1 

fl  =  --n(i-n),  f2  =  (i-rj)(i+n),  f3  =  -n(l+b)  (3.6) 

via 

N3J+I-3  =  fl(nl)  fjCfl2)  »  !»J  =  •  (3.7) 


The  coordinates  (x1?x2)  in  the  physical  space  and  the  displacements  (uru2)  are  given  by 
9  9 

xi  =  ^  NKXiK’  ui  =  ^  NKUiK  1  =  1’2  (3-8) 

K=1  K=1 


where  XiK  are  the  nodal  coordinates  and  are  the  nodal  displacements. 

T.  e  discrete  form  of  J  based  on  3x3  Gaussian  quadrature  appropriate  to  the  9-node 
biquadratic  Lagrange  element  is 


J  = 
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(3.9) 


The  field  quantities,  including  qx  and  its  spatial  derivatives,  are  evaluated  at  the  nine 
quadrature  points  and  weighted  by  Wp  and  the  determinant  of  Jacobian  matrix, 
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det  (3xk/3nk).  To  be  consistent  with  the  isoparametric  formulation  the  values  of  qx  and 
dqjdx^  at  the  quadrature  points  are  evaluated  from 

9  *  9  2 

3q  i  r  r  3Nj  3nk 

iT=  2  2  iTS7Ql1  '  (310) 

1=1  k=l 

where  Nj  is  the  biquadratic  shape  function  and  Qjj  is  the  value  of  qt  associated  with 
the  Ith  node  of  an  element.  Nodal  values  Qjj  are  assigned  in  accordance  with  a 
smooth  function.  Numerical  experiments  have  shown  that  the  value  of  J  is  insensitive 
to  the  type  of  smooth  function.  Thus  mesh  design  and  convenience  are  the  only 
considerations  in  selecting  a  smooth  function.  A  detailed  discussion  of  several  aspects 
of  the  implementation  of  the  domain  representation  of  J,  including  candidate  qj 
functions,  has  been  given  by  Shih  et  al.  (1986). 


3-5  Interaction  Integral  for  Extracting  Complex  Stress  Intensity  Factors 

Under  the  assumption  of  linear  elasticity  G  =  J,  and  therefore  G  has  the 
line-integral  representation 
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J°ikeik5lj 
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3x, 


n:  dr. 


(3.11) 


The  interaction  energy  release  rate  (see  Appendix  II,  eq.  (11.18))  is  defined  by 


^int  Gm  “  G 


(3.12) 


<7tot  is  the  energy  release  rate  of  the  total  field  (the  actual  field  plus  the  auxiliary  field) 
and  ^aux  is  the  energy  release  rate  of  the  auxiliary  field.  It  follows  from  (3.11)  and 
(3.12)  that 
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We  have  used  the  reciprocity  theorem  to  make  the  connection 


2  °ij(eij)aux  +  2  (ai.paux  eij  aij^€ij^aux  • 


(3.14) 


By  virtue  of  (3.12)  and  the  path-independence  of  the  integral  in  (3.11),  the  integral  in 
(3.13)  is  also  path-independent  for  an  elastic  homogeneous  medium  and  for  the  bimaterial 
problem  under  consideration. 

As  mentioned  in  Section  3.3,  it  is  more  advantageous  to  work  with  domain 
representations  in  finite  element  computations.  Using  the  weighting  function  qj 
discussed  in  Section  3.3,  the  domain  representation  of  (3.13)  is 
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(3.15) 


where  A  is  any  annular  domain  surrounding  the  crack  tip.  Replacing  the  terms  within 
the  brackets  [  ]  in  (3.9)  by  the  terms  within  the  brackets  [  ]  in  (3.15),  we  have  the  finite 
element  representation  for  G. 

int 

The  method  of  extracting  complex  stress  intensity  factors  will  be  explained  using  the 
particular  boundary  value  problem  under  study.  A  more  complete  discussion  is  given  in 
Appendix  II.  Let  (ojk)aux,  (€jk)aux  and  (Quj/Sxj)^  be  the  singular  fields  for  the  deformable 
medium  corresponding  to  unit  value  of  the  stress  intensity  factor  Kx  (designated  by  kt  in 
Appendix  II).  The  singular  stressses  are  given  in  Appendix  I.  The  strains  can  be  obtained 
using  li  t\C  s  law  and  du^/dx^  can  be  obtained,  after  some  effort,  from  available  solutions 
(e.g.,  Sih  and  Rice,  1964).  The  integral  in  (3.15)  is  evaluated  using  the  actual  numerical 
fields  oi j  and  Buj/axj,  evaluated  at  the  Gauss  points,  and  the  auxiliary  fields  (from  Appendix 
I)  also  evaluated  at  the  same  Gauss  points.  Let  Gj^  denote  the  value  of  the  interaction 
integral.  As  discussed  in  Appendix  II  (see  eqs.  (11.22)  and  (11.23)),  the  value  of  Kx  of  the  actual 
numerically  determined  fields  for  the  problem  at  hand  is  given  by 

2cosh^7re 
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Similarly,  with  the  help  of  auxiliary  singular  fields  appropriate  to  unit  value  of 
K2’  Glnt  is  evaluated.  The  value  of  K2  of  the  actual  numerical  fields  is  then  given  by 
2cosh^  77 € 


=  (2) 
^  (l-v)/fi  lnt 


(3.17) 


It  may  be  noted  that  the  method  can  also  be  applied  to  inclined  crack  and 
three-dimensional  crack  problems.  In  the  case  of  isotropic  homogeneous  bodies,  the 
usual  Kj,  Kjj  and  Km  stress  intensity  factors  are  extracted  from  the  numerical  fields  by 
using  auxiliary  fields  appropriate  to  the  isotropic  homogeneous  medium  and  the 
well-known  relation  between  the  energy  release  rate  and  Kj,  Kn  and  Km.  Such  an 
approach,  using  the  line-integral  (3.13)  to  calculate  what  is  in  effect  G.  „  has  been 
employed  by  Yau  et  al.  (1980),  and  Stern  et  al.  (1976)  for  isotropic  homogeneous  media. 


3.6  Computational  Model 

With  reference  to  the  geometry  depicted  in  Fig.  3a,  only  the  right  half  of 
deformable  medium  need  be  considered  in  the  finite  element  analysis  since  the  problem 
possesses  reflective  symmetry  with  respect  to  the  vertical  plane  bisecting  the  crack.  The 
half  crack  length  is  a,  and  the  half  width  and  height  of  the  deformable  slab  is  100a. 
The  finite  element  model  is  constructed  using  9-noded  quadrilateral  Lagrangian  elements 
(see  Fig.  3b)  .  We  employ  an  arrangement  of  wedge-shaped  9-node  elements  in  the 
immedi  crack  tip  region  -  short  of  embedding  the  actual  singularity  fields  into  the 
elements,  such  an  arrangement  of  elements  contains  sufficient  degrees  of  freedom  to 
reproduce  the  qualitative  features  of  crack  tip  fields. 

The  wedge-shaped  element  has  a  radial  length  of  10  ^a.  The  arrangement  of  wedge 
elements  is  surrounded  by  semi-circular  strips  of  elements  as  shown  in  Fig.  3c.  Each  decade 
of  radial  length  is  spanned  by  four  semi-circular  strips  of  elements  -  thus  the  domain  between 
10  ^a  and  a  is  spanned  by  60  strips  of  elements  which  are  generated  by  a  logarithmic  scale. 


Within  each  strip,  the  angular  distance  from  0  to  77  is  spanned  by  12  equally  sized  elements. 
The  mesh  for  the  domain  r  S  a,  as  laid  out  in  this  manner,  has  732  elements. 

The  domain  beyond  r  >  a,  bounded  by  the  remote  boundaries  and  the  symmetry  plane, 
is  modelled  by  140  elements.  All  in  all,  there  are  872  elements  and  3639  nodes  in  the  model. 
To  test  the  adequacy  of  the  mesh  for  the  problem  at  hand,  we  carried  out  several  elastic  and 
elastic-plastic  test  calculations  using  the  above  mesh  and  a  refined  mesh  which  had  twice  as 
many  elements.  The  results  differed  typically  by  less  than  1  percent  and  we  chose  to  use  the 
coarser  mesh  to  generate  the  results  to  be  reported  in  the  next  section.  The  adequacy  of  the 
mesh  design  and  the  element  type  in  conjunction  with  selective  integration,  and  the  accuracy 
of  the  domain  representation  for  calculating  J  and  the  interaction  integral  will  become 
evident  in  Section  4. 


4.  RESULTS 

4.1  Mixed  Mode  Crack  Tip  Fields  in  Homogeneous  Media 

The  plane  strain  elastic-plastic  fields  for  the  center  cracked  panel  have  been 
computed  for  the  case  where  the  upper  and  lower  half  materials  are  identical  and 
are  described  by  J2-deformation  theory  with  power  law  hardening.  In  this  case  the 
near  tip  fields  are  essentially  of  the  type  originally  derived  by  Hutchinson  (1968a, b) 
and  Rice  and  Rosengren  (1968)  for  pure  modes  I  and  II  and  by  Shih  (1974)  for 
mixed  n  ode.  Hereafter,  they  will  be  referred  to  as  HRR  fields.  Application  of  the 
J-integral  to  the  mixed  mode  small-scale  yielding  problem  reveals  that  the  dominant 
singularity  governing  the  asymptotic  behavior  of  the  stresses,  strains,  and 
displacements  near  the  crack  tip  has  the  form 
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In  the  above  expressions  the  dimensionless  angular  functions  o-,  and  uj 

depend  parametically  on  the  plastic  mixity  parameter,  Mp,  and  the  hardening 
exponent,  n.  The  plastic  mixity  is  defined  as 


Mp  =  -  tan"1 
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oee(r,e=0) 


/im 


r-0  ar0(r,9=O) 


(4.2) 


such  that  Mp  =  1  for  pure  mode  I  and  Mp  =  0  for  pure  mode  II.  The  amplitude  of  the 
HRR  singularity  field,  Kp  (plastic  stress  intensity  factor),  is  defined  by  Shih  (1974)  such 
that  the  angular  distribution,  oQ,  attains  a  maximum  value  of  unity.  With  this  definition, 
Kp  is  related  to  the  value  of  J  by 


aQo 

J  =  —  In(Kp)n+1  .  (4.3) 

The  factor  In  depends  on  the  degree  of  mixity,  Mp,  and  n,  and  has  been  determined 
for  a  wide  range  of  these  parameters  by  Shih  (1974).  For  our  present  purposes  we 
find  it  convenient  to  rescale  Kp  by  setting  In  =  1. 

For  strictly  linear  elasticity,  the  mixity  parameter  can  be  reinterpreted  as 


Mp  -  Me  ■  -  tan"1 
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For  small  scale  yielding,  where  the  stresses  beyond  the  plastic  zone  are  those  of  the 
elastic  K-field,  Me  is  also  defined  by 
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where  r*  is  within  the  zone  of  dominance  of  the  elastic  field.  In  this  case  J  is  given  by 
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For  the  case  of  a  finite  crack  in  an  infinite  homogeneous  plate,  two  types  of 
loading  are  analyzed,  one  for  which  the  elastic  K  field  corresponds  to  pure  mode  I, 
and  the  other  for  which  the  mode  is  mixed.  We  define  <p  as 


ip  =  tan 


<X) 

-ipi' 

a> 

^  n  * 
u22 


rK 


=  tan" 


II 


IKt 


=  77/2(1 -Me) 


(4.7) 


Thus  the  two  cases  correspond  to  ip  =  0°  and  -30°.  We  chose  to  display  the  results  for 
<p  =  -30°  because  the  mixed  plastic  field  that  develops  under  this  loading  bears 
noteworthy  similarities  to  the  field  near  the  interface  crack  tip  when  the  remote  loading 
on  the  bimaterial  plate  is  pure  tension. 

Figures  5a  and  5b  show  distributions  of  stress  ahead  of  the  crack  tip  in  a  nonlinear 
homogeneous  medium.  The  hoop  stress,  a6g,  and  effective  stress,  o&,  along  the  central  row 
of  quadrature  points  in  the  first  row  of  elements,  9  =  7.5  ,  are  shown  for  the  cases 
where  n  =  3  and  10.  The  radial  distances  are  normalized  by  rp,  the  length  of  the 
plastic  zone  along  0  =  0.  Note  that  the  stresses  plotted  are  the  actual  values  computed 


at  the  Gauss  points  -  we  do  not  apply  any  smoothing  to  the  data  points  in  this  and  other 
figures  in  this  paper. 

In  Fig.  5a  /n(stress/oQ)  is  plotted  on  the  ordinate  and  the  dotted  line  is  the  hoop 
stress  according  to  the  K  field.  The  regions  of  elastic  K-field  dominance  are 
indicated.  Within  the  plastic  zone  the  curves  indeed  take  on  the  slopes  -1/4  and  -1/11  in 


accorda  e  with  (4.1a).  In  Fig.  5b  the  normalized  stress,  stress/(a0[J/(aa0€0r)]^/( n+*) ), 
is  plotted  on  the  ordinate.  The  stress  levels  according  to  the  HRR  singularity  fields 


(4.1a)  for  n  =  3  and  10  are  indicated.  (These  levels  correspond  to  values  of  the 
dimensionless  stress  oqq  which  are  1.26  and  2.16  respectively.)  The  normalized  stresses  are 
nearly  invariant  for  r  <  rp  and  indicate  fields  that  are  entirely  of  pure  mode  I  HRR 
form.  Figures  6a  and  6c  show  the  angular  distribution  of  normalized  stress  (defined  as 
in  Fig.  5b)  deep  inside  the  plastic  zone,  for  n  =  3  and  10,  and  ip  =  0°.  Figures  6b 
and  6d  show  the  angular  variations  in  stress  for  the  mixed  mode  field  corresponding  to 
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i/>  =  -30°.  These  normalized  fields  which  are  obtained  from  a  full  field  analysis  agree 
perfectly  with  the  mixed  mode  asymptotic  fields  (Shih,  1974).  One  noteworthy  feature  of 
these  fields  is  that  the  levels  of  normal  stress,  for  a  given  overall  value  of  the  J-integral, 
are  reduced  for  the  mixed  mode  field.  This  is  not  only  the  case  for  the  normal  stress 
levels  directly  ahead  of  the  crack,  but  is  also  true  for  the  maximum  values  of  stress  over 
the  full  range  of  9.  Another  interesting  feature  of  the  mixed  mode  solutions  is  the  large 
shear  stresses  ahead  of  the  crack  tip. 

Plastic  zones  determined  from  the  analysis  of  the  full  boundary  value  problem 
are  shown  for  these  same  values  of  n  in  Figs.  7a  and  7b  for  pure  mode  I  loading 
(i/»  =  0°)  and  mixed  mode  loading  (0  =  -30°).  The  plastic  zones  are  essentially 
identical  to  those  obtained  by  Shih  (1974)  based  on  an  analysis  of  the  small  scale 
yielding  boundary  value  problem.  For  the  mixed  mode  case,  the  plastic  zones  are 
no  longer  symmetrical  above  and  below  the  crack  line  and,  for  a  fixed  value  of 
remote  loading  as  measured  by  the  value  of  J,  extend  further  than  for  the  pure 
mode  I  case  (note  the  difference  in  the  scale  of  the  plots). 

Contours  of  hydrostatic  stress  are  shown  in  Figs.  7c  and  7d.  Here  again,  it  is 
noteworthy  that  with  pure  mode  I  loading  the  contours  are,  as  expected,  symmetric  about 
the  crack  line  but  become  highly  skewed  under  mixed  mode  loading.  A  feature  of 
the  hydrostatic  stress  field  for  the  n=10  material  is  worth  attention:  at  a  fixed  overall 
load  level  as  measured  by  the  value  of  J,  the  hydrostatic  stress  contour,  =  oQ,  for  the 
mixed  mode  field  extends  further  than  the  contour  for  the  pure  mode  I  field;  however, 
the  zone  of  high  hydrostatic  stress,  oh  2aQ,  is  smaller  under  mixed  mode  field. 

4.2  The  Interface  Crack  in  Linear  Elastic  Media 

Figure  8  shows  results  for  the  near  tip  stress  fields  of  the  interface  crack  for 
v  =  0.3.  In  Figs.  8a-8d  normalized  stresses  are  plotted  in  the  form  stress/(JE/r)^.  The 
analytically  obtained  asymptotic  field  (see  Appendix  I)  is  shown  by  the  dotted  line  and 
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the  essentially  exact  correspondence  between  it  and  our  numerical  solution  is  evident.  At 
distances  r/a  <  0.05  from  the  crack  tip,  the  differences  between  the  numerical  fields  of 
the  full  boundary  value  problem  of  the  center  cracked  plate  and  the  analytic  asymptotic 
fields  are  less  than  5  percent  (see  Fig.  8a).  This  establishes  a  zone  of  dominance  for 
the  asymptotic  linear  elastic  field. 

Figure  8a,  like  Fig.  5,  shows  the  radial  distribution  of  stress  along  a  ray 
connnecting  the  central  row  of  quadrature  points  in  the  elements  along  the  bond  line. 
Note  that  the  field  is  of  a  strongly  mixed  nature  characterized  by  large  shear  stresses. 
As  noted  in  Section  2.1,  the  stresses  begin  to  oscillate  as  the  crack  tip  is  approached  and, 
oqq  first  becomes  compressive  at  r/a  =  10"9  at  9  =  7.5°. 

Figures  8b-8d  show  angular  variations  of  stress  at  three  radial  distances  from  the 
crack  tip.  The  strong  dependence  of  the  angular  variations  of  the  fields  on  the 
relative  radial  distance  is  evident. 

To  extract  the  stress  intensity  factors  Kx  and  K2  from  the  numerical  fields,  the 
interaction  integral  (3.15)  was  evaluated  on  several  semi-circular  strips  of  elements  with 
mean  radii  ranging  from  10'14a  to  a.  The  values  of  <7jnt  obtained  from  the  various 
strips  differed  only  in  the  third  or  fourth  significant  digit.  Specifically,  for  v  ranging 
from  0.0  to  0.499,  the  value  for  and  K2,  as  determined  by  (3.16)  and  (3.17),  agreed 
with  the  analytical  solution  given  by  Rice  and  Sih  (1965)  (see  Fig.  2)  to  better  than  2 
significant  digits.  The  excellent  agreement  attests  to  the  quality  of  the  numerical  fields 
and  the  accuracy  of  the  domain  interaction  integral  (3.15)  for  uncoupling  and  evaluating 
stress  intensity  factors. 

4-3  Small  Scale  Yielding  for  the  Interface  Crack  Subject  to  Remote  Tension 

As  discussed  previously  a  characteristic  feature  of  the  asymptotic  linear  elastic 
solution  is  the  increasingly  rapid  oscillations  of  the  fields  as  the  crack  tip  is  approached.  The 
relative  crack  face  displacements,  given  by  (2.5),  imply  wrinkling  of  the  faces  and 
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interpenetration  near  the  tip  of  the  crack.  The  finite  element  solutions  accurately 
reproduce  the  oscillation  in  the  stress  field  (see  Fig.  8)  and  crack  face  penetration  into 
the  rigid  substrate  for  distances  greater  than  10‘15a  from  the  crack  tip. 

However,  using  the  nonlinear  model  based  on  a  power  law  hardening  behavior,  we 
find  that  at  remote  load  levels  which  are  a  fraction  of  the  yield  stress,  the  crack  face 
opens  up  smoothly.  The  finite  element  computations  gave  no  indications  that  penetration 
into  the  rigid  substrate  would  develop  even  at  distances  smaller  than  10'15a.  The  stresses 
increase  monotonically  as  we  probe  deeper  into  the  plastic  zone  and  there  is  no  pattern 
that  suggests  that  the  stresses  would  reach  a  peak  and  then  decrease  as  the  distance 
within  the  plastic  zone  becomes  vanishingly  small.  For  one  particular  load  level  we 
repeated  the  entire  calculation  using  a  mesh  which  had  twice  as  many  elements.  The 
numerical  fields  agreed  to  within  two  significant  digits,  strongly  suggesting  that  fields  at 
distances  of  10‘10a  to  10'15a  are  accurately  resolved  by  our  computations. 

At  each  stress  level  the  value  of  J  was  extracted  from  the  numerical  fields  using 
(3.9)  and  the  'plateau'  qj  function.  The  values  of  J  evaluated  on  annular  domains  with 
mean  radii  ranging  from  10'14a  to  a  generally  agreed  to  within  3  significant  digits.  At 
stress  levels  where  the  maximum  plastic  zone  size  is  much  smaller  than  the  crack  length, 
the  value  of  J  agreed  precisely  with  the  analytical  result  of  (2.16),  viz.. 


J  =  G  = 


(K*  +  K2). 


(4.8) 


4  cosh27le 

The  essentially  exact  path-independence  of  the  numerically  evaluated  J  and  the  precise 
agreement  with  (4.8)  under  small  scale  yielding  demonstrate  the  quality  of  the  numerical 


solution  and  the  accuracy  of  the  domain  method.  The  J  values  thus  obtained  are 
employed  for  the  purpose  of  normalizing  the  fields  to  be  discussed. 

We  choose  to  present  results  for  a  high  hardening  material  (n=3)  and  a  moderate  to 
low  hardening  material  (n=10)  to  reveal  the  effect  of  plasticity  on  the  near  tip  fields. 
In  particular,  the  behavior  of  the  hoop  stress  ahead  of  the  crack  is  shown  in  detail  since 
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it  is  rather  descriptive  of  the  conditions  at  the  interface  and  also  plays  an  important  role 
in  fracture  processes.  Also,  the  variation  of  the  hoop  stress  with  radial  distance  is 
representative  of  the  other  stress  components.  In  particular  the  hoop  stress  normalized 
by  the  remotely  applied  stress  provides  a  direct  indication  of  the  intensification  of  the 
stress  field.  As  before,  the  hoop  stress  at  points  on  a  radial  line  passing  through  the 
central  row  of  Gauss  points  in  the  elements  bonded  to  the  interface  (9  =  7.5°)  is  chosen 
for  detailed  examination.  Since  the  distances  under  discussion  range  from  10~15a  to  a  it 
is  necessary  to  plot  log(r/a)  on  the  absicca.  Over  distances  that  differ  by  15  orders  of 
magnitude,  the  normalized  hoop  stress  varies  by  several  orders  of  magnitude,  and  in  some 
cases  is  compressive.  With  the  exception  of  a  small  interval  where  the  hoop  stress  goes 
to  zero  and  changes  sign,  the  magnitude  of  the  normalized  hoop  stress  is  much  larger 
than  unity  within  the  distances  under  discussion.  To  grasp  an  overall  picture  of  the 
behavior  of  the  fields,  we  confine  our  attention  to  normalized  stresses  with  magnitudes 
greater  than  unity.  This  permits  the  results  to  be  presented  in  a  rather  compact  and 
informative  way.  Specifically  we  take  the  log  of  |a00/a°°|  and  to  preserve  the 
algebraic  sign  of  the  stress,  the  negative  sign  is  appended.  Accordingly,  we  plot 
sign(a0e)log  |  o00/ o°°|. 

The  normalized  stress  determined  by  the  finite  element  calculations  for  the  strictly 
elastic  material  is  plotted  against  normalized  distance  in  Fig.  9a  (see  dashed-line  curve). 
The  break  in  the  dashed-line  curve  near  log(r/a)  =  -9.0  corresponds  to  the  stress 
changing  from  tension  to  compression.  It  should  be  noted  that  the  position  of  the  break 
in  the  dashed-line  curve  is  independent  of  the  remote  tensile  stress  level.  We  now 
examine  the  behavior  of  the  stress  for  the  strongly  hardening  (n=3)  nonlinear  material. 
At  the  lowest  load  level  considered,  o“/oQ  =  2.0xl0~5,  the  plastic  zone  is  confined  to  a 
distance  of  about  10"na.  While  the  stress  is  negative  within  the  plastic  zone,  the  slope  of 
the  curve  becomes  positive  for  r/a  <  10'14.  At  the  next  load  level,  o°°/a0  =  2.0xl0'4,  the 
stress  is  positive  over  the  entire  distance  studied,  however,  at  the  outer  fringe  of  the 
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plastic  zone,  10  9  <  r/a  <  10'8,  the  trend  of  the  surrounding  elastic  field  is  felt.  At 
distances  well  inside  the  plastic  zone  the  hoop  stress  increases  monotonically  as  the  crack 
tip  is  approached.  At  a  still  higher  load  level,  o  / aQ  =  6xl0~3,  the  hoop  stress  increases 
monotonically  over  the  entire  distance  under  discussion  —  there  is  no  trace  of  an 
oscillatory  field. 

The  behavior  of  the  hoop  stress  for  a  moderate  to  low  hardening  material  is  shown 
in  Fig.  9b.  For  comparison  purposes,  the  linear  elastic  result  is  again  included.  At  the 
very  low  remote  stress  level  a°°/o0  =  2.0xl0‘5,  the  stress  in  the  outer  plastic  zone  is 
compressive  in  response,  we  note,  to  the  compressive  pressure  of  the  surrounding  elastic 
fields.  Well  within  the  plastic  zone  the  stress  is  tensile.  At  the  next  remote  stress 
level,  a°°/o0  =  2.0xl0'4,  the  hoop  stress  in  the  outermost  fringe  of  the  plastic  zone  is 
weakly  affected  by  the  elastic  fields.  A  short  distance  into  the  plastic  zone,  the  hoop 
stress  increases  monotonically  as  the  crack  tip  is  approached.  At  the  final  load  level, 
°  /°o  =  6.0x10  3,  shown  in  the  plot,  there  is  no  trace  of  any  oscillatory  field.  (The  solid 
line  is  terminated  at  log(r/a)  =  -12  because  we  were  unable  to  complete  the  calculations 
for  the  n=10  material  at  this  load  level  using  a  mesh  which  could  resolve  fields  at 
distances  of  10  a  from  the  crack  tip.)  We  also  note  that,  even  at  the  largest  remote 
stress  level,  the  nonlinear  crack  tip  region  is  surrounded  by  an  annular  zone, 
10'5  <  r/a  <  10'2,  in  which  the  stresses  are  well  approximated  by  the  singular  fields 
(2.13). 

To  develop  a  better  understanding  of  the  plastic  fields  we  detail  their  angular 
distribution  at  two  radial  distances  within  the  plastic  zone  and  at  three  load  levels: 
a“/ao  =  2x1  O'4,  om/a0  =  6xl0'3  and  a°/aQ  =  2X10'1.  We  begin  with  the  n=3  material.  The 
length  of  the  plastic  zone  as  measured  along  the  interface  is  denoted  by  rp.  The  relative 
plastic  zone  size,  rp/a,  corresponding  to  the  three  remote  loads  are  5.6xl0'8,  3.2xl0‘5  and 
l.OxlO'2  respectively.  As  mentioned  previously,  the  asymptotic  structure  of  the  nonlinear 
crack  tip  fields  has  yet  to  be  identified  -  thus  we  do  not  have  a  definitive  form  with 
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which  to  organize  the  numerically  calculated  fields.  For  reasons  which  will  be  evident 
shortly,  the  stresses  are  normalized  using  an  HRR  structure;  in  particular,  the  normalized 
stresses  are  explicitly  defined  as: 


A 


°ij/ao 

(J/aa0e0r)1/(n+1) 


(4.9) 


The  angular  distribution  of  the  normalized  stresses  at  a  fixed  relative  distance 
r/rp  =  lxlO'2,  are  shown  in  Figs.  10a,  c  and  e  for  the  three  remote  load  levels.  The 
angular  fields  deep  inside  the  plastic  zone,  r/rp  =  l.OxlO-6,  are  shown  in  Figs.  10b,  d 
and  e.  It  is  clear  from  these  plots  that  the  fields  are  quite  similar  to  the  mixed 
mode  HRR  fields,  for  positive  values  of  0,  shown  in  Fig.  6b.  Specifically,  within 
any  given  plastic  zone  the  stress  fields  shift  toward  a  mode  I  HRR  type  angular 
distribution  as  the  distance  moves  deeper  into  the  plastic  zone  (e.g.  compare 
Figs.  lOe  and  f).  At  a  fixed  distance  relative  to  the  plastic  zone  size,  the  angular 
fields  also  shift  toward  a  mode  I  HRR  type  angular  distribution  as  the  remote  load 
or  the  size  of  the  plastic  zone  increases  (e.g.  compare  Figs.  10c  and  e).  Finally  we 
point  out  that  within  the  range  of  distances  and  loads  examined,  the  normalized 
stresses  are  of  order  unity;  the  maximum  value  of  the  normalized  hoop  stress  vary 
between  0.8  and  1.4. 

It  is  instructive  to  examine  the  changes  in  the  size  and  shape  of  the  plastic  zone  as 
the  remote  load  increases.  To  facilitate  comparisons  with  the  usual  presentation  of  these 
plots  (e.g.  Shih,  1974),  the  distances  are  normalized  by  JE/{(1 -v2)o2}.  The  plastic  zones  at 
the  three  load  levels  are  plotted  in  Fig.  11a,  using  these  normalized  distances 


X  =  x(1-v2)o2/JE  and  Y  =  y(l -v2)a2/JE  . 


(4.10) 


The  relative  size  of  the  plastic  zone  decreases  (when  plotted  on  these  non-dimensionalized 
axes)  and  changes  over  to  a  nearly  mode  I  pattern  as  the  remote  load  increases. 
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Hydrostatic  stress  contours  for  o ^  -  a0  and  =  2.0aQ  are  shown  in  Figs,  lib,  c  and  d. 
The  zone  of  high  hydrostatic  stress,  £  2o0,  increases  in  size  as  the  field  shifts 
towards  a  mode  I  distribution. 

We  detail  the  fields  for  the  n=10  material  at  remote  stress  levels,  o*/oQ  of  2xiCT4, 
6x10  3  and  2X10'1.  The  respective  plastic  zone  sizes,  rp/a,  are  8.2xl0'8,  4.2xl0~5  and 
1.2xl0"2.  Plots  of  the  angular  fields  at  distances  of  r/rp=  lxlO'2  and  lxlO'4  for  the  three 
remote  stress  levels  are  shown  in  Fig.  12.  It  is  easily  seen  that  the  angular  distribution 
of  the  plastic  fields  are  nearly  those  of  the  mixed  mode  HRR  type  and  that  the  fields 
deep  inside  the  plastic  zone  (as  in  Fig.  12f)  shift  toward  the  mode  I  HRR  distribution  as 
shown  in  Fig.  6c. 

Plastic  zones  corresponding  to  the  three  remote  stress  levels  are  shown  in  Fig.  13a. 
Again  the  shift  towards  a  mode  I  pattern  with  increasing  remote  stress  is  easily  seen. 
Hydrostatic  stress  contours  are  plotted  in  Figs.  13b,  c  and  d;  note  the  increase  in  the  size 
of  the  zone  of  high  hydrostatic  stress,  ah  £  1.7o0,  in  Fig.  13d. 

5.  DISCUSSION 

Our  results  suggest  that  a  finite  element  method  which  uses  selective  integration 
(Malkus  and  Hughes,  1978;  Hughes,  1980)  and  9-node  Lagrangian  quadrilateral  element  is 
well  suited  for  studying  the  complex  nonlinear  stress-strain  fields  near  the  tip  of  an 
interface  crack.  We  find  the  method  to  be  highly  stable;  the  reduced  integrations  of 
volumetric  fields  using  4  quadrature  points  display  no  hour  glass  or  other  spurious  modes. 
The  method  is  easily  programmed  and,  we  note,  lends  itself  to  coding  which  is  readily 
vectorized  for  supercomputers. 

The  full  field  numerical  results  presented  in  Section  4.3  indeed  suggest  that  the 
near  tip  fields  of  the  interface  crack  do  not  have  a  separable  form  of  the  type  given  by 
the  HRR  singularity  (4.1).  Nevertheless  our  small  scale  yielding  solutions  for  the 
interface  crack  bear  remarkable  similarities  to  the  fields  for  cracks  in  homogeneous 


-30- 


media  subject  to  mixed  mode  loading.  Furthermore  within  the  distances  that  we  have 
examined,  the  elastic-plastic  fields  do  not  exhibit  any  of  the  non-physical  incompatible 
displacements  and  oscillatory  stresses  that  are  characteristic  of  the  small  strain  linear 
elastic  solutions.  We  discuss  these  features  in  turn. 

For  the  boundary  value  problem  at  hand  we  note  that,  under  small  scale  yielding 
conditions, 

(a2”)277a(l+4e2) 

J  =  - Z - (l-v)/<t.  (5.1) 

4coshz  Tie 

If  JQ  is  the  value  of  J  for  a  similar  crack  in  a  homogeneous  plate  then,  for  a  fixed 
remote  stress  a^2, 

j  (l+4e2) 

y  = : — •  (5.2) 

o  2cosh  77€ 

For  the  case  of  the  crack  on  a  rigid  substrate  with  v  =  0.3  in  the  deformable  medium, 
e  =  0.0935  as  already  noted,  and  J/JQ  =  0.475.  We  now  recall  the  distributions  in  each 
of  Figs.  6,  10  and  12  where  the  stresses  are  normalized  by  the  quantity 
Q  =  a0[J/(°£ao€or)]1  ^n+1')  ■  With  this  perspective  it  is  clear  that  for  a  fixed  intensity  of 
remote  loading,  as  measured  by  o"2,  the  stresses  near  the  tip  are  generally  lower  for  the 
interface  crack  than  those  for  the  crack  in  a  homogeneous  media.  However,  at  a  given 
value  of  J  the  stresses  near  the  tip  of  the  interface  crack  are  generally  higher;  at  larger 
distances  from  the  crack  tip,  the  stresses  are  comparable  or  lower  than  those  for  the 
homogeneous  media.  For  example,  the  maximum  value  of  (ogg/q)  in  a  homogeneous 
medium,  under  mixed  mode  loading  corresponding  to  ip  =  —30°,  is  ogg/q  =  1.1  when  n=3 
and  oqq/ q  =  1.7  when  n=10.  For  the  interface  crack,  the  maximum  stress  levels,  at 
each  radial  distance,  depend  on  distance  and  remote  stress  level.  At  the  normalized 
distance  r/rp  =  l.OxlO'2,  Ogg/q  =  0.9,  1.1,  and  1.25  for  the  three  remote  stress  levels 
plotted  in  Fig.  10  n=3.  At  the  distance  r/rp  =  l.OxlO'6,  ogg/q  =  1.2,  1.35,  and  1.45 
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which  are  considerably  larger  than  those  for  a  crack  in  a  homogeneous  medium.  The 
same  trends  hold  for  n=10.  It  should  be  noted  that  the  angular  distribution  of  the 
stresses  for  the  interface  crack  can  be  quite  different  from  those  for  the  crack  in  a 
homogeneous  medium.  An  important  difference  is  the  large  shear  stresses  and  strains 
that  develop  near  the  bond  line  of  an  interface  crack  despite  the  tensile  nature  of  the 
remote  stresses. 

Another  intriguing  feature  of  the  stress  fields  of  the  interface  crack  is  how  the 
degree  of  mixity  of  the  near  tip  fields  changes  with  the  extent  of  plasticity.  For 
example,  comparing  Figs.  12  and  6c  and  6d  for  n=10,  it  may  be  seen  that  as  the  remote 
stress  level  increases,  and  as  the  positions  examined  are  closer  to  the  crack  tip,  the  field 
more  nearly  resembles  that  of  a  pure  mode  I  HRR  field.  Indeed,  near  tip  fields  for 
remote  stress  o*/oo  >  0.2  are  essentially  of  pure  mode  I  distribution.  A  discussion  of 
these  fields  is  deferred  to  a  later  paper.  At  lower  stress  levels,  and  at  distances  more 
remote  from  the  crack  tip,  the  field  has  features  similar  to  the  homogeneous  mixed  mode, 
for  example  that  for  ip  =  “30°. 

As  noted  in  Sections  2  and  4,  the  small  strain  linear  elastic  asymptotic  field  for  an 
interface  crack  is  characterized  by  displacement  incompatibilities  and  oscillatory  stresses. 
In  Section  2  it  was  noted  that  under  combined  remote  tension  and  shear,  the  regions 
ahead  of  the  crack  where  the  stresses  oscillate  can  be  significant  fractions  of  the  crack 
length,  thus  rendering  the  solutions  physically  unacceptable.  Within  the  nonlinear  theory 
of  plane  stress  Knowles  and  Sternberg  (1983)  have  carried  out  an  elegant  asymptotic 
analysis  of  an  interface  crack  between  two  neo-Hookean  sheets.  Their  asymptotic  fields, 
based  on  nonlinear  kinematics  and  a  linear  relation  between  the  Cauchy  shear  stress  and 
shear  displacement,  are  free  of  oscillatory  singularities;  they  also  found  that  the  crack 
opens  smoothly  near  its  ends.  Our  full  field  solutions  show  that,  even  within  the 
framework  of  linearized  kinematics,  accounting  for  nonlinear  material  behavior 
effectively  precludes  incompatible  displacements  and  stress  oscillations.  This  is  true  not 
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°nly  for  the  case  of  remotely  applied  pure  tension  but,  as  we  have  noted,  for  the  more 
general  case  of  combined  remote  tension  and  shear.  To  demonstrate  this  we  analyzed 
the  case  where  a22  >  0  and  o12  =  “2o22,  or  as  described  in  connection  with  (2.11)  and 
(2.12)  where  X  =  -58.7°.  The  results  shown  in  Figs.  14  and  15  were  obtained  by 
performing  a  boundary  layer,  small  scale  yielding  analysis  in  which  the  stresses,  given  by 
the  asymptotic  field  in  Appendix  I,  were  imposed  at  a  distance  which  is  large  compared 
to  the  plastic  zone  size.  The  stress  intensities  for  this  case  were  prescribed  as 
=  °2 2^  i'/tta  and  K2  =  o2 2{2+2e}v,77a  with  o22  taken  to  be  the  scalar  stress 
variable  (see  Fig.  2c). 

Figure  14a  shows  the  linear  elastic  results  for  the  stress  variations  along  a  radial 
line  at  0  =  7.5  for  the  combined  remote  tension  and  shear  loading.  It  may  be  recalled 
from  Section  2  that  on  the  bond  line  oqq  is  compressive  at  r/a  =  2xl0"3.  At  0  =  7.5°, 
Oq9  is  compressive  at  r/a  =  3x10  4.  Figure  14b  shows  the  radial  distributions  of  stress, 
again  along  the  central  row  of  quadrature  points  at  0  =  7.5°,  for  a  nonlinear  material 
with  n=3  and  for  three  remote  stress  levels.  In  the  plot  the  hoop  stress  oQQ  is 

normalized  by  o22  which  is  labelled  as  o°°.  Note  that  with  combined  loading  of  this  sort, 
that  is  with  cr”2  <  0,  larger  remote  levels  of  o“2,  as  compared  to  the  cases  dominated  by 
remote  tension,  are  required  to  'open  the  crack'  and  produce  monotonically  increasing 
stresses  as  the  crack  tip  is  approached.  Nevertheless  even  at  the  rather  low  remote 
stress  lc  cl  of  <^22/  <"'o  =  1-0X10  4  oscillatory  stresses  and  overlapping  of  the  crack  faces 
are  precluded  in  a  high  hardening  material.  In  a  moderate  to  low  hardening  material, 
the  oscillatory  fields  are  precluded  at  substantially  lower  remote  tensile  stresses.  Indeed, 
even  at  rather  low  stress  levels,  the  characteristics  of  the  plastic  fields  (which  are  totally 
different  from  those  predicted  by  linear  elasticity)  prevail  over  significant  length  scales. 

Figure  15  shows  the  angular  distribution  of  stresses  at  six  radial  distances  within 
the  plastic  zone  for  the  remote  loading  corresponding  to  o22/oo  =  l.OxKT1  and 
aiJao  =  -2-OxlO'1.  For  this  load  level,  rp/a  =  5.0xl(T2  which  is  also  the  maximum 
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extent  of  the  plastic  zone.  A  remarkable  feature  of  these  fields  is  the  gradual  shift 
towards  a  mode  I  like  HRR  field  as  the  distance  r/rp  becomes  vanishingly  small  despite 
the  relatively  large  remote  shear  stresses.  This  trend  can  be  clearly  seen  by  comparing 
the  angular  variations  of  the  fields  in  Fig.  15  to  the  fields  in  Fig.  10  which  are 
associated  with  pure  remote  tension. 

Despite  the  pathological  nature  of  the  linear  elastic  singular  fields,  we  must  point 
out  that  there  is  a  finite  annular  region  where  the  elastic  singular  fields  (2.13)  is  a  good 
approximation  of  the  full  field  solution  as  can  be  seen  in  Fig.  9  and  Fig.  14.  It  is  in 
this  more  restricted  sense,  that  we  speak  of  small  scale  yielding;  m  the  present  context  the 
term  does  not  imply  self-similar  growth  of  the  plastic  zone.  In  a  sequel  to  this  paper,  the 
near  tip  fields  for  an  isolated  crack  and  collinear  arrays  of  cracks  under  contained 
yielding  and  essentially  fully  yielded  conditions  will  be  presented  and  implications  for 
interface  fracture  will  be  discussed. 
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APPENDIX  I 


Singular  Fields  Near  the  Tip  of  an  Interface  Crack 

As  noted  in  the  text,  the  asymptotic  crack  tip  fields  for  the  case  of  linear  isotropic 
media  take  the  form 


a=  1 


(a) 

”^ZZ  gij  [0,6 /n(r/2a);e]  . 
V2m 


(1.1) 


The  stress  intensity  factors,  Ka,  in  (1.1)  were  defined  in  (2.1)  with  examples  given  in 

Fig.  2.  For  the  bimaterial  geometries  shown  in  Fig.  2  the  full  form  of  the  singular 

(oc) 

field  in  region  1  is  given  by  (1.1)  with  the  gjj  defined  as  follows.  Let  (ij)  denote  the 
polar  indices  (r9)  and  define  A  as 


Then 


A  = 


(1) 

grr 


1 

2  cosh(7ie)  ’ 

=  A  |^e  €(77'®)  j*3  cos(0/2  +  e/n(r/2a))  +2e  sin9  cos(0/2  -  e/n(r/2a)) 

-  sin0  sin(0/2  -  e/n(r/2a))| 

-e^77'0)  cos(39/2  +  €/n(r/2a))l  ,  (1.3) 


(2) 

grr 


(1) 

§00 


A  ^e'6^77"0)  |-3sin(0/2+  e/n(r/2a))  +  2e  sin0  sin(0/2  -  e/n(r/2a)) 
+  sin@  cos(9/2  -  e/n(r/2a))j 
+  e€(77'0)sin(30/2  +  e/n(r/2a)) 

A  e  €(77'®)  jcos(9/2  +  e/n(r/2a))  -2esin0  cos(0/2  -  e/n(r/2a)) 


(1.4) 


+  sin0  sin(0/2  -  e/n(r/2a)) 


+  e^77’0) 


cos(30/2  +  e/n(r/2a)) 


(1.5) 
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gee  =  A  |e"e(n‘0)  |-  sin(0/2+e/n(r/2a))  -  2e  sin9  sin(0/2  -  e/n(r/2a)) 

-  sin9  cos(9/2  -  €/n(r/2a))j 

-  ee(n~  0)  sin(30/2  +  €/n(r/2a))J  ,  (1.6) 

=  A  |^e‘€(7T‘9)  jsin(0/2  +  €/n(r/2a))  -2e  sin0  sin(0/2  -€/n(r/2a)) 


-sin9  cos(0/2  -e/n(r/2a))j 

+  ee(n~e)  sin(30/2  +e/n(r/2a))j  ,  (1.7) 

s[2Q  =  A  [c"^77-0)  jcos(0/2  +  e/n(r/2a))  +  2e  sin9  cos(0/2  -e/n(r/2a)) 

-sin0  sin(0/2  -e/n(r/2a))  j 

+  e6^71"0)  cos(30/2  +  e/n  (r/2a))j  .  (1.8) 


These  asymptotic  fields  can  easily  be  rearranged  into  the  form  given  by  Sih  and  Rice 
(1964)  by  redefining  K  as  their  k  as  described  in  Section  2. 
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APPENDIX  II 


Crack  Field  Mode  Interaction 


In  this  section  we  describe  a  method  for  calculating  the  individual  stress  modes  for 
a  crack  subjected  to  mixed  mode  loading.  The  method  is  applicable  to  cracks  in 
isotropic  or  anisotropic  media.  Our  method  is  motivated  by  Eshelby’s  (1956,  1961)  notion 
of  interaction  forces  and,  in  the  present  context,  this  means  an  interaction  energy  release 
rate,  C7int.  To  begin,  consider  the  field  of  a  finite  size  crack  (such  as  shown  in  Fig.  1)  in 
an  arbitrary  anisotropic  but  homogeneous  elastic  media.  The  remote  tension  loading,  a"2, 
is  augmented  with  the  remote  shear  stress  a”2  (as  depicted  in  Fig.  2c).  Let  xr  x2  be  the 
Cartesian  coordinates  such  that  the  line  crack  lies  along  the  xx-axis  over  the  region 
~a  ^  xx  $  a.  Along  the  crack  line 


ai2  ai2 


(x2-  a2) 
i 


J4  ’ 


|xi  I  >  a,  x2  =  0  ,  i  =  1,  2,  3  . 


(II.  1) 


At  a  vanishingly  small  distance  from  the  right  crack  tip  and,  with  5  defined  as  in  the 
text  (Fig.  1), 

a 

,  5  >  a,  x2  =  0  .  (II. 2) 


ai2  -  °i2 


25 


Note  that  on  the  crack  line  the  tractions  are  independent  of  elastic  constants  and  are 
essentially  set  by  static  equilibrium  (Barnett  and  Asaro,  1972). 

In  homogeneous  media  the  notion  of  stress  intensity  factors  is  tied  to  the  symmetry 
in  the  stress  fields  which  are  expressed  in  the  form  (r,0  are  polar  coordinates  centered  on 
the  crack  tip) 


III 


°ij  = 


u  (a) 

Ka/(27lr)*f>:'  (0) 


a=I 


(II.3) 


where 


-37- 

4(9=0)  =  1  ,  4(0)=  o  =  4(0)  (II.4  !> 

4(0=0)  =  1  ,  4(°)  =  0  =  4(°)  (IL42) 

and 

hi  hi  hi 

f32(0=O)=l  ,  f22(0)  =  0  =  f12(0)  (I1.43) 

For  the  problem  at  hand  then, 

oo  - 

K i  =  a22  w  »  KII  =  a-  ✓an  ,  Kra  =  o”  ✓a n  .  (II-5l,2,3) 

The  displacement  fields  will  not  possess  such  simple  symmetries,  and  again  for  the 
problem  at  hand  we  note  that  the  jump  in  displacement  across  the  crack  tip  is  given  by 
Barnett  and  Asaro  (1972)  as 

Aui  =  4  aj2  ✓a/n  vx/2n  ,  |x|  S  a  (II.6) 

where  B  ^  is  the  inverse  of  the  so-called  prelogarithmic  energy  factor  matrix,  B. 

For  plane  strain  we  can  define  a  stress  intensity  vector  with  components 

K  =  {Kn  ,  Kj)  (H.7) 

such  that  the  crack  traction  vector  and  the  displacement  jumps  on  the  crack  line 

1  =  {°L2’  a22>  (II.8) 

and 

Au  =  {AUp  Au2)  (H.9) 

are  given  by 

1 

t  =  -  K  ,  x  >  a  (II.  10) 

v2nx 

and 

1  -1  - 

=  —  B  •  K  -/x/277  ,  |  x  |  <  a  . 


Au 


(11.11) 
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The  energy  release  rate  can  be  calculated  from  the  virtual  work  integral  of  (2.14) 
and  is  given  by  Barnett  and  Asaro  (1972)  as 

1  -1 

C  =  J^KocBcc  0KB  (11.12) 

where  the  indices  ot,B  can  be  interpreted  such  that  1  is  identified  with  II  and  2  with  I. 
Note  that  for  isotropic  media  B  is  diagonal,  with  components 


“  B22  ~ 


- - —  and  B„„ 

4  77(  1  -v)  33 


477  ’ 


(11.13) 


and  the  usual  expressions  for  t,  Au,  and  G  are  recovered,  i.e.. 


G  = 


(1-v) 

2k 


r  2  2  n  1 

[K.  +  K..] +  yy  K 


2 

III 


(11.14) 


Stress  mode  decoupling  in  a  given  mixed  mode  stress  field  is  accomplished  through 
computation  of  what  might  be  referred  to  as  an  interaction  energy  release  rate.  To 
appreciate  the  procedure,  the  details  of  calculating  the  Kx  component  of  a  total  crack  tip 
field  are  illustrated. 

We  first  note  that  for  the  actual  (mixed  mode)  crack  field  the  crack  extension  force 
is  given  by 


c  ■  IT  [  K‘  +  2KiS’i*K*  + 


13  ^3  + 


(terms  not  involving  Kx) 


(11.15) 


Now  suppose  that  an  auxiliary  pure  mode  I  field,  of  intensity  k.,  is  superimposed  on 
the  actual  field  and  the  crack  extension  force  again  computed;  the  result  is 

(Kx  +  kx)  5ji( K1  +  kx)  +  2(K,  +  kj)  B\\  Kx 

+  2(K1  +  kx)  B1Z  K3  +  (terms  not  involving  or  kx) 


(11.16) 
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The  auxiliary  field  itself  has  the  energy  release  rate 

^  1  -1 
Caux  “  87j  kl  ’ 

From  the  above  relations  we  can  define  an  interaction  extension  force  as 


(11.17) 


Cint  =  Gtot  -  G  -  Caux=  ^2kj  B\ !  Kj  +  2k,  b\\  K1+  2kx  B*  Ks]. 

(11.18) 

All  three  energy  release  rates,  G,  Gtot,  and  Gaux,  including  Cjnt  itself,  are  computable  from 
either  line  or  domain  integrals  as  explained  in  Sections  3.3  to  3.5.  When  the  above 
calculations  are  performed  for  the  other  two  modes  a  simple  linear  system  of  equations  of  the 


form 


(a)  1  -l 

Cint  =  ~^T  ka  BasKs  ’  (no  sum  on  a  =  1,  2,  3) 


(H.19) 


results  which  may  easily  be  solved  for  the  Ks;  when  the  ka  are  assigned  unit  values  the 
solution  for  the  K$  are 


„(s) 


Ka  477  B Qg  Gint 


(11.20) 


It  may  be  noted  that  Cint  is  readily  evaluated  with  high  precision  via  the  domain 
integrals  described  in  Section  3.5. 

A  completely  analogous  procedure  can  be  applied  to  evaluate  the  respective 
intensities  of  the  two  stress  modes  K1  and  K2  for  the  present  interface  crack  being 
investigated.  For  plane  strain  we  note  by  (2.16)  that,  if  G  is  written  in  the  form 

-l 

G  =  Ka  Bcx8  K0  »  (11.21) 

that  B  is  diagonal  with  the  nonzero  components 


{(l-v1)//x1)/(4cosh277e) 


(11.22) 


Once  the  two 


are  evaluated,  the  two  stress  intensities  are  determined  from 


K 


a 


1 

2 


B  G 


(«) 

int 


(11.23) 
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FIGURE  CAPTIONS 


Figure  1: 
Figure  2: 
Figure  3: 

Figure  4: 
Figure  5: 

Figure  6: 

Figure  7: 

Figure  8: 

Figure  9: 

Figure  10: 

Figure  11: 

Figure  12: 

Figure  13: 


Interface  crack  between  two  bonded  dissimilar  materials. 

Stress  intensity  factors  for  three  crack  geometries. 

(a)  Interface  crack  between  deformable  medium  and  rigid  substrate. 

(b)  Finite  element  mesh  of  right  half  of  deformable  medium. 

(c)  Arrangement  of  elements  at  crack  tip. 

Conventions  for  domain  J.  Domain  A  is  enclosed  by  r,  C+,  CL  and  C 
Unit  normal  mj  »  nj  on  CQ  and  mj  =  -nj  on  r. 

(a)  Stresses  ahead  of  crack  tip  in  homogeneous  medium*  the  zone  of  K. 
dominance  is  indicated.  (b)  Normalized  stress  versus  distance 
normalized  by  plastic  zone  size. 

Angular  variations  of  normalized  stresses  in  homogeneous  medium  subject  to 
pure  mode  I  and  mixed  mode  loadings;  (a)  and  (b)  n=3  material;  (c)  and  (d) 
n=10  material. 

(a)  and  (b)  Plastic  zones  in  homogeneous  medium  subject  to  mode  I  and 
mixed  mode  loadings;  (c)  and  (d)  Hydrostatic  stress  contours  corresponding 
to  mode  I  and  mixed  mode  loadings. 

(a)  Normalized  stresses  (from  linear  elastic  calculations,  n=l)  near  the  bond 
line  versus  log  distance;  the  asymptotic  stresses  are  included,  (b),  (c)  and 

(d)  Angular  variations  of  normalized  stresses  at  three  radial  distances. 

Plots  of  log  of  normalized  hoop  stress  (near  the  bond  line)  against  log  of 
normalized  distance  for  the  interface  crack;  (a)  n=3  material,  (b)  n=10 
material.  The  elastic  (n=l)  result  is  included. 

Plots  of  angular  variations  of  normalized  stresses  for  the  interface  crack 
with  n=3  at  two  distances  within  the  plastic  zone  and  for  three  remote 
stress  levels. 

(a)  Plastic  zones  for  the  interface  crack  with  n=3  for  three  remote  stress 
levels,  (b),  (c)  and  (d)  Hydrostatic  stress  contours  for  three  remote  stress 
Levels.  The  normalized  distances  are  defined  as  X  =  x(l-v2)o2/JE  and 

Y  =  y(l-v2)o2/JE. 

Plots  of  angular  variations  of  normalized  stresses  for  the  interface  crack 
with  n=10  at  two  distances  within  the  plastic  zone  and  for  three  remote 
stress  levels. 

(a)  Plastic  zones  for  the  interface  crack  with  n=10  for  three  remote  stress 
levels,  (b),  (c)  and  (d)  Hydrostatic  stress  contours  for  three  remote  stress 
l_evels.  The  normalized  distances  are  defined  as  X  =  x(l-v2)o2/JE  and 

Y  =  y(l-v2)a2/JE. 
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Figure  14:  Hoop  stress  near  the  bond  line  due  to  combined  remote  tension  and  shear. 

(a)  Normalized  stresses  from  linear  elastic  calculations  (n=l);  the  asymptotic 
stresses  are  included,  (b)  Normalized  hoop  stress  for  nonlinear  material 
behavior  with  n=3.  The  elastic  (n=l)  result  is  included. 

Figure  15:  Plots  of  angular  variations  of  normalized  stresses  for  the  interface  crack 
with  n=3  at  six  distances  within  the  plastic  zone  for  remote  stresses  given  by 

°2  2/ao  =  °-1  and  or*/*. =  _o-2- 
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